Nonuniqueness and Structural Stability of Self-Consistent Models of Elliptical Galaxies

In recent years, high-resolution ground-based and HST observations have shown that the space density profiles ρ(r) near the central regions of most elliptical galaxies exhibit power-law cusps, i.e. they scale as ρ ∼ r with 0 ≤ γ ≤ 2. This finding contradicts the previously widely held view that ellipticals have extended central “cores” of near-constant density. Furthermore, evidence has accumulated over the past two decades that the constant-density surfaces of at least some ellipticals may actually be genuinely triaxial ellipsoids rather than spheroids. This has spurred some interest in the galactic community because theoretical and numerical work suggests that a central density cusp or mass concentration embedded in a generic nonaxisymmetric potential can lead to a significant expansion of the phase-space regions where motion is chaotic. In view of this evidence, Merritt and Fridman (1996) used Schwarzschild’s (1979) method to construct self-consistent models of a triaxial generalization of Dehnen’s spherical potential, which contains a central density cusp. They found that self-consistency could be achieved only for models with a weak (γ = 1) cusp, and then only when the chaotic orbits populating the outermost regions of the model were allowed to be not fully mixed (i.e., they were not required to sample the invariant measure [Lichtenberg & Lieberman, 1992]). Based on this and later work, Merritt and collaborators have suggested that triaxiality may not in fact be compatible with central density cusps, and that, as a consequence, most elliptical galaxies may be (or may be evolving towards) axisymmetric configurations, at least in the central regions. The construction of a self-consistent equilibrium ρ(r) via Schwarzschild’s numerical method is effected by assigning appropriate weights to a large number of orbital templates, each evolved under the influence of the gravitational potential generated by ρ(r), so that the weighted superposition of all the templates reproduces, at some level of coarse graining, the initial ρ(r). This is usually implemented via some constrained optimization algorithm, such as linear or quadratic programming. Owing to its conceptual simplicity and relative ease of implementation, Schwarzschild’s method has been used quite extensively for the construction of stellar equilibria, especially when no known analytical solutions to the collisionless Boltzmann equation exist or they are difficult to compute. In the present work (see also Siopis, 1999; Siopis & Kandrup, 1999) Schwarzschild’s method was first used to construct self-consistent models of a Plummer sphere and to calculate a number of velocity moments, which were then compared with known analytical solutions (Dejonghe, 1986) to assess the reliability of the numerical method. Subsequently, the method was applied to the construc-